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This is a pedagogical review of the lattice study of finite density QCD, which is intended 
to provide the minimum necessary contents, so that the paper may be used as the first 
reading for a newcomer to the field and also for those working in nonlattice communities. 
After a brief introduction to argue why finite density QCD can be a new attractive subject, 
we describe fundamental formulae which are necessary for the following sections. 

Then we survey lattice QCD simulations in small chemical potential regions, where 
several prominent works have been reported recently. Next, two-color QCD calculations are 
discussed, where we have a chance to glance at many new features of finite density QCD, 
and indeed recent simulations indicated quark pair condensation and the in- medium effect. 
Tables of SU(3) and SU(2) lattice simulations at finite baryon density are given. 

In the next section, we make a survey of several related works which may be a starting 
point of the new development in the future, although some works do not attract much 
attention now. Materials are described in a pedagogical manner. Starting from a simple 2-d 
model, we briefly discuss a lattice analysis of NJL model. We describe a non-perturbative 
analytical approach, i.e., strong coupling approximation method and some results. Canonical 
ensemble approach instead of usual grand canonical ensample may be another route to reach 
high density. We examine the density of state method and show this old idea includes 
recently proposed factorization method. An alternative method, complex Langevin equation 
and an interesting model, finite isospin model, are also discussed. 

In the Appendix, we give several technical points which are useful in practical calcula- 
tions. 
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§1. Introduction 



Finite density QCD (Quantum Chromodynamics) has attracted considerable 
attention in high energy physics, nuclear physics and astrophysics. Many theorists 
now believe that QCD has a very rich structure when we study it in temperature 
and density parameter space, and some experimentalists want to reveal it. Fig^ 
shows a schematical phase diagram in (T, fi) plane. Lattice QCD has been expected 
to provide fundamental information about QCD in nonperturbative regions, and 
recently many promising works have appeared in this field. 

In this paper, we discuss these activities together with a pedagogic introduction 
and related works. This is not a textbook-type comprehensive review. Minimum 
necessary knowledge and possible hints for future are given. Materials are taken 
from our notes, which we have been preparing during our research in this field. We 
hope this paper is a useful starting point for those who wish to join the field. For 
further reading, we cite several reviews .^''^''^'EJ''^ 

The most attractive possibility of finite density QCD is probably the color super 
conductivity (CSC). This was revealed first by one gluon exchange type calculation: 
Barrois, Bailin and Love, and Iwasaki and Iwado pointed out that an attractive force 
near the Fermi surface creates a Cooper pair, and results in CSC in the case of QCD 
at very low temperature and finite density.'^ In the late nineties, using the instanton 
type modeling of the attractive force, Alford, Rajagopal and Wilczek^' and Rapp, 
Schaefer, Shuryak and Velkovskjl^ argued that the gap energy is of the order of 
100 MeV, and therefore the transition temperature, which should be around the gap 
energy, is relatively high. In these arguments, two-fiavor, i.e., u and d quark, color 
superconductivity is considered and termed 2CS. Figj^jis a prediction of the phase 
by NJL model.l^ 

Since these works, the field has been actively explored and there have been 
many interesting observations, e.g., at sufficiently large chemical potential, a state 
with a special combination of three flavors and color is stable and called color flavor 
locking (CFL)^ a new phase LOFF (Larkin, Ovchinnikov, Fulde and Ferrell), where 
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pairing with different cliemical potentials results in nonzero total momentum, may 
be favored;'^^ ' the gap energy is ~ exp(— c/5) instead of ~ exp(— c/^^l^ and so on. 
See Ref. I13|) . Introduction of Ref. I14() also gives a good overview of the field. 




Fig. 1. Schematic phase diagram of QCD in {T,jj,) plane. 
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Fig. 2. Phase diagram {T, /i) plane by NJL model. (Ref. O). 

Another interesting, but probably less well known possibility is chiral restoration 
in the nuclear medium. In the normal nuclear matter, the baryonic density is po ~ 
0.17/fm^. In a Nambu-Jona-Lashino (NJL) type model, the quark condensation, (qq) 
at pq is less than that of the vacuum. We show the behavior of {i/^ip) as a function 
of /X and T in a simple NJL model calculation in Fig|31 The quark condensation 
is an important parameter and causes many effects, such as vector meson masses. 
There are several recent experimental data which may yield information regarding 
in-medium hadrons, 

• CERES: Large low mass e~^e~ pair enhancement was observed in CERN SPS in 
Pb+Au collisions at 158 A GeV/c. This non trivial enhancement is also found at 
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40 AGeV/c where it is even larger, and the data may only be reproduced if the 
properties of the intermediate p in the hot and dense medium are modified.^ 

• KEK: At KEK, invariant mass spectra of electron-positron pairs were measured 
in the region below the uj meson mass for the p+C and p+Cu collisions. The 
result was interpreted as being a signature of the p/u modification at a normal 
nuclear-matter density.^ 

• STAR: The invariant mass of p meson decays in Star experiment at ^/s = 200 
GeV Au+Au collisions at RHIC shows 60-70 MeV downward shift of the peak 
from its vacuum value. This suggests the modification of the spectral function 
at finite T and p^ 

• Pionic Is states of Sn: Precision spectroscopy of pionic Is states of Sn nuclei 
suggests /*(/>o)^//^ = 0.65 ± 0.05, at the normal nuclear density.^ 

See Ref. il9|) for more detailed arguments and related theoretical and experimental 
works. 



1 




Fig. 3. {qq) as a function of {T,fi). We solve gap equations at finite temperature and density 
derived in Ref. I2UII . We employ all parameter values given in I2UII . except tGo, which is set to 
be zero, i.e., we drop the s-quark contribution for simplicity. The arrow indicates the value of 
jj, corresponding to the normal nuclear density po- 

Hagedorn noteci^ that if the hadron spectrum increases exponentially, 

/9(m) = Ce'"/^^, (M) 

then the partition function diverges at some temperature Tc, called Hagedorn tem- 
perature, 

J dmp{m)e-"'/'^ = oo if T > Tq. (1-2) 
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Cabbibo and Parisi pointed out that the singularity of the partition function does 
not necessary mean the upper hmit of the temperature, but the phase transition!^ 

In 1981, McLerran and Svetitsky,!^ and Kuti, Polonyi and Szlachanyi,!^ demon- 
strated that when the temperature increases, a QCD matter encounters a phase tran- 
sition from the confinement to the deconfinement phase by calculating Polyakov lines 
for two-color QCD. Since these pioneering works, lattice simulations have provided 
many detailed analyses of QCD at finite temperature. Still the order of the phase 
transition, i.e., the cross-over, the first or the second order, with the physical it, d 
and s quark masses is not yet finally determined. Lattice QCD is the first principle 
calculation, and can be considered as the most reliable starting point. See Ref. 05)1 
for recent lattice investigations at finite temperature. 

Then how about lattice QCD at finite density ? A lattice simulation of color 
SU(2) at finite temperature and density was carried out in 1984 and the transi- 
tion from the confinement to the deconfinement state was observed.!^ Since then, 
to our knowledge, few lattice calculations were performed until 199^^'!^ except 
algorithmic study. 

The reason is that when the chemical potential is introduced, the fermion de- 
terminant, det /\(/u) becomes complex, which appears in the Euclidean path integral 
measure, 

Z = Tle-^^^-^'^^ = j VUVi}V^e-^^'^-^^'^ = j VU dei Ae'^^^ . (1-3) 

Here A is diagonal in flavor, A = diag(Z\u, Z\^, Ag, ■ ■ ■ )■ 

The simplest way to avoid the problem is to use the quench approximation, where 
we discard det A{^). This approximation works well in many lattice calculations such 
as spectroscopy. However at the finite density it was found that the onset of chiral 
symmetry restoration occurs at a chemical potential of half the pion mass, 

/^c = ^ (1-4) 

i.e., at the chiral limit, /ic = 0.^^' Fig|l]shows data points in the quark mass-chemical 
potential plane at which physical observables start deviating from their = values 
in quench simulations at /? = 0. The curve (dotted line) indicates the expected 
critical line correpoinding to pion (baryon) threshold. 

Stephanov pointed out that quenched QCD is not the A'^^ limit of QCD, 
but the limit of a theory with 2Nf quarks.!^ Let us consider the resolvent of the 
Dirac operator A for zero mass case, 

G{x,y) = (Tr^-) = / dx' dy' p{x\y')^-^, (1-5) 
z — A J z — z' 

where z = x + is a complex variable. Banks-Casher relation is ('0V') = "^/"(O)- P 
can be inversely obtained as 

P=-^- (1-6) 



vr 



dz* 
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Fig. 4. Phase dil];ram in the quark mass-chemicll potential plane of the SU(3) ^luge theory with 
staggered fermions at /3 = in the quench approximation!^ The data points show the threshold 
value where physical observables start deviation from their fi — values. The dotted line the 
expected critical line corresponding to the baryon threshold and the curve is the pion threshold 
which seems to describe the Monte Carlo data. The dotted curve is the result of the mean field 
calculation. 



G 



Using the formula, det A = exp Tr log A, we can write G as 

dV 
dz ' 
where 

= -(logdet(z-Z\)). 
If the standard replica trick works, V can be obtained as lim„^o where 

1 

n 



Vn = --\og{det^{z- A)). 
Taking tI)l,r = (1 i ^^)/2'tJj as bases, we can write A as 

A = 







iX + fi 




(1-7) 
(1-8) 

(1-9) 

(MO) 



where we employ the chiral representation of the Dirac matrices. Using the random 
matrix model where X is treated as complex Gaussian random variables, Stephanov 
found that the above Vn does not lead to the natural result, and we should start 
from Vn which includes not only z but also z* , 



n 



log(det"(z-Z\)(z* -Z\t)). 



Ml 
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This corresponds to a system composed of n quarks with original action and n 
quarks with the conjugate action. Therefore to understand QCD at finite density, 
the quenched approximation is not appropriate and one needs dynamical simulations. 
We have learned : 

• QCD at finite density is an interesting attractive field of physics, and lattice 
QCD should play an important role here. 

• Quench approximation of lattice QCD is not appropriate for the finite density 
system. 

• Due to the complex nature of the fermion determinant detZ\, the standard 
Monte Carlo simulation is very difficult to obtain physical quantities. 

In the following we continue our discussion, assuming the reader accepts these points. 

In Sec. 2, we give the fundamental formulae in lattice QCD at finite baryon 
density. In Sec. 3, we describe recent prominent activities in small chemical potential 
and high temperature regions. In Sec. 4, several numerical simulations of two-color 
QCD are reported, which suggest new interesting features at large /U. In Sec. 5, we 
collect several ideas and calculations which may give us a hint for the next step. We 
summarize and conclude our discussions in Sec. 6. 



§2. Formulation 



A thermodynamical system is described by the partition function, 

Z = Tr(e-w(-f^-M^)) 

where /5 = 1/T and we impose the periodic (anti-periodic) boundary condition for U 
(V') along the temporal directiorF^ . The chemical potential is introduced in fermion 
action as ^tp^Qip. Therefore the fermion determinant has the form, 

Z\(/i) = D^^^ + m + ^70. (2-2) 

Adopting anti-hermite D and hermite 7 (final result does not depend on the 
representation.), we can easily show 

Z\(^)^ = -Dy-^y + m + /i7o = 75Z\(-/i)75 (2-3) 

and then 

(det A{^)y = det A{^i)^ = det Z\(-^). (2-4) 

At /i = 0, this relationship guarantees that det A is real. 

In the continuum field theory, the chemical potential is introduced by the sub- 
stitution, 

P4 ^ P4 - i/^ (2-5) 

in the fermion propagators. On the lattice with the lattice spacing a, the ipd^ip term 
leads to 'il>e'^^'^°"ip. For example, a free propagator of Wilson fermions is given as, 

1 



1 - ^ Ej=i {(1 - 7^.)e'P^'' + (1 + 7/.)e-*P^'^} ' 



(2-6) 
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Therefore, we can naturally include the chemical potential exponentiated in the 
fermion matrix as, 

3 

A{x,x') = 5,,,, - K {(1 - 70f^.(^)5.',.+J + (1 + 7i)Uj{x')6^,^,_i} 

i=l 

_^ |e+M(l - 74)C/4(x)5,, + e-f^il + 74)^1 .4} • (2-7) 

For staggered (Kogut-Susskind) fermions, 

1 ^ 

Aix,x') = m5,y + - ^r?,(x) {u.{x)5^,^^^i - 

i=l 

+^r?4(x) {e+'^C/4(x)5,, ,^4 - e-^?7i(x')5,,,,_4} , (2-8) 

where 

r^.ix) = 1 = 1), (-1)-^ (/. = 2), (-1)-^+-^ = 3), (-l)-i+-^+-3 = 4). 

(2-9) 

In other words, the chemical potential is introduced as, 

Utix) ^ e^Utix), 

Uhx)^e'^u}{x), (2-10) 

where U^{x) = exp(iaA^(x)). Hasenfratz and Karsch have shown that the formula 
(|2-lUp avoids nonphysical divergence of the free energy of quarks Gavai considered 
more general form than exp(ib/i) in Ea. ()2-10p .-^ Creutz discussed how the chemical 
potential appears in lattice fermion formulationP^ 

The relation H2-3|) holds for Wilson fermions, while in staggered fermion case, it 
is easy to check 

z\(/x)t = ^^A{-^Ji)^^, (2-ii) 

where 

n{x, x') = (-l)-i+-2+x3+x45^^^,. (2-12) 

In the relativistic formulation, we expect the invariance under — > —ji. If we 
expand the partition function as 



^(^) = E^|^ / ^f/detZ\e-/^^-U=o, (2-13) 



j^k Qk 

< I 

k 

the terms with odd power of ^ vanish. Let us check k = 1 case explicitly; At ^ = 0, 



d dA dA 

— detA = Tr[A-^—] det A = TV[75Z\-S575^75] det75Z\75 

= TY[Z\t-i(-M^)]detZ\t = -Tr[A{U^)-^ ^^^] det A{U^), (2-14) 
o/j, dfi 
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Fig. 5. Eigen value distribution. Wilson fermions. 4^ lattice, ma — 0.1. Quench (3 — 5.7. 

where we use Ea. ()2-3() . The integration measure VU and the gluon action Sq are 
invariant under the change ?7 — > f7^, we get dZ/dfi = —dZ/d^ = at /i = 0. 

We show in Figs El and IHl typical eigen value distributions for Wilson fermions 
(Fig[2I) and staggered fermions (FigE)). Here configurations are generated in the 
quench approximation, since so far there is no algorithm which enables us to calculate 
the path integral (|l-3() at finite /i. 

It is sometimes convenient to change variables, V and ip as 

ip{x,ti) e"*'^^(f,ti) 

i>{x,ti) ^ e+'^^4>{x,t,) (2-15) 

where = 1, 2, • • • , Nt. Then 

ti)e^Ut{x, ti)ip{x, ti+i) 4;{x, ti)Ut{x, ti)ip{x, ti+i) (2-16) 

i.e., we do not need to put a factor e^^ to the link variables, Ut and uj , except at 
the temporal edge of the lattice where 

e^^''Ut{x,Nt) = e''/^Ut{x,Nt) 
e-^'f'Ujix, Nt) = e^f'/^Ujix, Nt). (2-17) 

The partition function should be a function of Ntfj. = = n/T (see Ea. (|l-3|) ). 
therefore the above expressions appear to be natural. 
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§3. Three Color 

Even though the fermion determinant det A is complex, one may perform Monte 
Carlo simulations by taking its modulus as a measure, 

-03g 

(O) = ^ I VUO det A e-'^^'^ '"^ 



1 

Z 



f VUO\detA\e^^ e- 
J PC/| detZ\|e*^ e" 



(3Sg 



fVUO\ det A\e'^ e'^^G / VU\ det A\e'^ e 



id p-pSa 



jVU\detA\e-P^G ' jVU\detA\e 



(3-1) 



where (• • • )o is the expectation value with | det A\ as the measure, i.e., phase quench- 
ing measure. The direct calculations were pursued on small lattices, but a large phase 
fluctuation hinders us from obtaining a meaningful signal in low temperature and 
large chemical potential regions, i.e., the numerator and denominator of the last 
term of Ea. (|3-1() are very small.'^''^''^ Lattice QCD at finite chemical potential 
suffers from a severe sign problem. In spite of this difficulty, there have been many 
challenging efforts, which we will survey in this section. Table 01 is a compilation of 
numerical simulations of three color system. 

3.1. Response of observables with respect to fi at fj, = 

Although the direct simulations at finite /i is very hard, one may measure the 
effect of the chemical potential through the response of physical observables with 
respect to the chemical potential at /i = 0. Such an attempt was first pursued by 
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Action 


parameters 


comments 


Ref. 


Plaq. + KS 


Nf = 2, m, = 0.00625,0.0125, 
16^ X 8 


quark number susceptibility 




Plaq. + KS 


Nf = 2, = 0.0125, 0.0170, 0.0250, 
16 X 8^ X 4 


response of observables at = 




Plaq. + KS 


Nf = 2+1, m„,d = 0.025, = 0.2, 
N^ x4, N,^ 8,10,12 


phase diagram, q method 




Imp. + Imp. KS 


Nf = 2, m, = 0.1,0.2, 16^ x 4 


Taylor expansion, reweighting 
method 


EJ 


Imp. + Imp. KS 


Nf = 2, m/T = 0.1, 0.4, 16^ x 4 


Taylor expansion, analytic frame 
work, equation of state 




Plaq. + KS 


Nf = 4, = 0.05, 16^ x 4 


imaginary chemical potential 




Plaq. + KS 


Nf = 2, m, = 0.025, 8^ x 4, 6^ x 4 


QCD phase diagram, imaginary 
chemical potential 




Plaq. + KS 


Nf = 2+1, mu,d = 0.025, = 0.2, 
iVf x4,Ns = 8,10,12 


equation of state 





Gottlieb et. al."^ 

They calculated the singlet and non-singlet quark number susceptibilities, Xsi+) 
and XNs{-), 

XS,NS = 7^ — ± ^ — [riu ± rid) (3-2) 



where for i = u,d is the quark number densities given by 

TdlogZ , , 

ni = — ^ (3-3) 

V dm 

These susceptibilities are found to be very small in the low temperature phase and to 
increase abruptly at the transition temperature. The two susceptibilities are almost 
same, e.g. xs ~ Xns which indicates that fluctuations in the u and d quark densi- 
ties are uncorrelated. These susceptibilities are important since they are related to 
event-by-event fluctuations and diffusion in relativistic heavy-ion collisions.!^ Fur- 
ther studies in the susceptibilities for quenched and 2 flavors QCD can be found 
[j^MMiJ^EB MILC (MIMD Lattice Computation) collaboration recently calcu- 
lated the susceptibilities for three flavors with improved staggered fermions.^^' FiglT] 
shows the triplet susceptibility Xtrip 

(= XNS here) on 8^ x 4, 12^ x 6 and 16^ x 8 lat- 
tices and indicates the excellent scaling properties between the Nt = 6 and Nf = 8 
results. FigElshows the difference between singlet and triplet susceptibility. One can 
see the clear difference between the two susceptibilities around the phase transition 
which occurs at about 180 MeV. 

QCD-TARO (Thousand-cell ARray processors for Omni-purposes) collaboration 
has extended the idea of Ref.'^^' further and calculated derivatives of meson screening 
masses and the quark condensation with respect to the chemical potential at = 
Consider a meson propagator C{z) in z direction consisting of u and d quarks as an 
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100 



400 



200 300 
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Fig. 7. The triplet quark number susceptibility fo Nf = 3 on 8^ x 4, 12^ x 6 and 16^ x 8 lattices. 



example, 



C{z) = A{e-^' + e 



Mz _|_ M{L^-z)\ 



(3-4) 



= ^ y VUH{z)H^O) e-^-^G-^^^ 

= ^ j VUG{z) deiAe-^^^, 

where H{x) = '4){x)r'il){x) is a meson operator and 

G{z) ^T,rA-\z,{))rA-/{Q,z). 

We neglect possible higher states for simplicity. 

The first derivative of C{z) with respect to the chemical potential leads to 



(3-5) 



(3-6) 



1 dCiz)_ldA^m^(^_L.^^^^^ 



C{z) dfj, A dfi dfi 
Similarly the second derivative is given as 



1 d^C{z) _ 1 d^A f 2 dAdM ^ d^M' ^ ^ 



M z 



C{z) d^j? A dfi"^ \Adfj, dn dfi^ 



z 1 tanh 

2 



h 

2 



(3-7) 



2 
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Fig. 8. The difference between singlet and triplet quark number susceptibility for Nf = 3 on 8 x 4, 
12^ X 6 and 16^ x 8 lattices. 



The numerical data of the derivative of C{z) to is given as follows. 

dC{z) d < G{z) > ' 



dn 



(3-9) 



< Gjz) > 
dp? 



-2{G + G 




■(G) 




(3-10) 

where D stands for det A and the dotted O and O denote the first and the second 
derivatives of an operator O with respect to p. Fitting the numerical data of Eas. H3-9() 
and (|3-l(jp to Eas. (|3-7j ) and (|3-8|) we can determine the derivatives of the residue A 
and the meson mass M. Thus with the derivatives of M, we can obtain the behavior 
of a hadron mass in the vicinity of /i = as a function of the chemical potential, 



M{^i) = M(0) + iJL 



dM 



dp 



/i=0 



1 ^d^M 



+ 



(3-11) 



The first derivative of the pseudo-scalar (PS) mass turned out to be consistent 
with zero ^ = 0. In Fig|Hl we show the second derivative of the PS mass where 
iso-scalar and iso-vector type chemical potentials, ps and ^y, are introduced. 



IJ-V = fJ'u = —fJ'd- 



(3-12) 
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In the low temperature phase the second derivative of the PS mass M with respect to 
fis is small. This is expected since below the critical temperature and in the vicinity 
of zero fiSi the PS meson still keeps a property as a Goldstone boson and persists 
its behavior. On the other hand, above Tc, the PS meson is no longer a Goldstone 
boson and the second derivative of M seems to remain finite. In contrast to the case 
of iJ,s the second derivative of M with respect to fiy is significantly large in the low 
temperature phase and decreases in magnitude above Tc. 
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Fig. 9. The second derivative of the pseudo scalar mass for isoscalar and isovector chem- 
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ical potentials, /is and nv. (Ref. I39ll l 



As we will see later, the fermion determinant with the isovector chemical poten- 
tial is equivalent to that of a two-flavor system without the phase. Large nonzero 
values of the second derivative of the pseudo-scalar meson mass may suggest some 
phase structure in finite but small chemical potential regions. Son and Stephanov 
have shed a new light on the model as a finite isospin density system. See section 

Gavai and Gupta calculated the depedence of the pressure P on /x through a 
Taylor expansion at fj, = For a homogeneous system, the pressure P is obtained 
by P = —F/V, where F = — Tlog Z. Therefore the Taylor expansion of P at /ij = 
{i = u,d,...) is given by 

1 (9F 1 d^F 
P(,„. ...) = P(0) - _ ^ _ + (3.13) 

i id 

The second term Xi corresponds to the quark number density which is zero at /i = 
and the third therm Xiji the quark number susceptibility. One can consider the 
higher order derivative terms and all the odd derivatives vanish by CP symmetry. 
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The Taylor series of P up to the fourth order for = = ^ is given by 

2' 



AP 

rp4 



where AP = P{ii) - P(0), and 



J^2 



(3-15) 



(3-16) 



Xuudd is neglected since it turns out to be numerically smaller than the errors in 
Xuuuu- FigOHshows obtained in quenched QCD. For /i/Tc = 0.15, the effects of 
/X on P are very small, which implies that at RHIC region where ///Tc ~ 0.06, the 
effects are also small. As n/Tc increases, the effects become significant. One could 
see the effects at SPS region where ///Tc ~ 0.45. 




AP 

Fig. 10. as a function of T/Tc for the values of ji/Tc. 



3.2. Glasgow approach 

The Glasgow group has been developing the way towards simulations of lattice 
QCD at finite density. See Ref. I54|) and references therein. 

Consider the partition function at /i normalized by that of /i = 0. 

jVUdetA{^,)e-P'o / det Z^(0) e-/^^^ 

jVUdetA{{))e-P^G j VU det A{id) e-P^c 



Taking the derivative of Z with respect to T and we obtain the energy and the 
baryon number density, respectively. By investigating zeros in the complex ^-plane, 
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Z{fi) = (3-18) 

we have Lee- Yang zero. It was very difficult to obtain reliable values of {b\n\) numer- 
ically at low temperature when n increases. 

3.3. Reweighting by Fodor and Katz 

Fodor and KatJ^ have achieved a great step in lattice QCD simulation at finite 
density, i.e., they have for the first time succeeded in estimating the phase transition 
line at finite /U in (T, n) parameter plane. See Figllll 

They proposed using the multi-parameter reweighting method. The essential 
idea is to use the gauge coupling constant /? as a controlling parameter, 



(O) 



-i- [ PC/0 detZ\(^)e-^^G 

1 . /■wO^^^4^e-(^-^°)^«detZ\(0)e-*^« 
J detZ\(0) ^ ^ 



,dcl 
dct zi(0) 



\^ dot A(0) ^ I 



\ det zi(0) ^ I 

and choose 3n so that 



det Z\(/u) 



detZ\(0) 



(3-19) 



(3-20) 



becomes as large as possible, i.e., the overlap should be large. In Ea. ()3-19() , (• • • )o 
stands for the evaluation at {Po,fJ, = 0). The coupling constant /?o is taken at the 
confinement/deconfinement transition point, and physical quantities near the phase 
transition line are calculated by the formula (|3-19() . i.e., quantities on the phase 
transition line in (T, plane are evaluated at (Tc, = 0) with the reweighting factor. 
A possible interpretation for this procedure is that points on the phase transition 
line have similar nature. fFig |12() 

If the reweighting factor, (|3-2Up . is large, we can get better signal. Further if 
in some regions of {fi, (3) parameter space, this value changes little, we have similar 
reweighting effect there. Ejiri investigated the problem, and observed that the chem- 
ical potential effect through the phase of the fermion matrix does not contribute to 
the (3 dependence of the gauge action, and in the amplitude, | det /^(//)|, correlates 
strongly with e"'^^'^^ ,S6i) 

In this calculation we need the ratio of the determinant, det Z\(/i)/Z\(0). Gibbs 
obtained a useful formula,*^ 

detZ\(^) = e^'^=^JJ(Ai + e-^*^), (3-21) 
1=1 



where Aj are eigen values of a matrix constructed from the spatial and mass terms 
of fermion matrix, which does not depend on the chemical potential ^. Therefore 
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Fig. 12. Schematic diagram of multi-parameter reweighting. 



once we obtain the eigenvalues A,, we can calculate the determinant at any /i for a 
fixed gauge field. See Appendix A. 

Note that, though Fodor and Katz obtained great success near the phase tran- 
sition line, they did not solve the sign problem of the finite density lattice QCD. It 
still remains difficult to study a low-temperature region where the phase fluctuation 
coming from the complex fermion determinant is most probably large. Fig |13l shows 
(cos 9) as a function of /x obtained by Bielefeld-Swansea group for various beta and 
quark masses. At large ^, the complex phase of the fermion determinant becomes 
large, i.e., the nominator and denominator in Eo 13- II are small. 

3.4. Taylor expansion 

Fodor and Katz calculated determinants in Eq. (|3-19j) . but it is difficult to study 
large lattices by the present computers. In any case, we cannot study large chemical 
potentials because of the sign problem. 
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Fig. 13. The phase of complex determinant calculated by Bielefeld-Swansea collaboration. 
(Ref-Ell). 



Bielefeld-Swansee collaboratiod^ proposed using Taylor expansion at /i = 0, 
f detAii,) \ _^/." 9"lndetZ\ 

\detA{0)J ~ ^^nl ^'^ > 

At the lowest order, the phase of the determinant is obtained as ulmTrM"^— — 

a/i 

which is easier to evaluate than the direct calculation from the determinant itself*^. 

Employing the formula, they have performed thorough studies in small ji regions 
which cover RHIC experiments and obtained the phase diagram similar to that of 
Fodor and Katz. 

Recently they study the thermodynamical grand potential up to the fourth order 
of the derivative with respect to the chemical potential, and calculate the quation of 
state. |42") Fig llSI shows A{p/T^) which is defined as 

(3-23) 



\ _ p 




v_ 








rp4 


T,0 



The correction of the pressure is large for 0.9 < T/Tq < 1.3, h/Tq > .5, but will 
decrease as T rises further. Fodor, Katz and Szabo reported similar result. I43|l 



TrM — — is usually estimated by the noise method which uses noise vectors but this term 
might need a number of vectors to obtain its value with accuracjl^ 
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Fig. 14. Phase diagram obtained by Bielefeld- Swansea collaboration. (Ref. I4lll '). 

3.5. Imaginary chemical potential 

If the chemical potential is pure imaginary, i.e., ^ = z^/, then the fermion matrix 
A becomes, 

and 



A = Du^y + m + i;U/7o 



(3-24) 
(3-25) 



Z\t = 75Z\75. 

Thus det A with a pure imaginary chemical potential is real. 

If we use the expressions 12-171 in section 2, the link variable at t = Nt has the 
form, 

e'^Ut{x,Nt). (3-26) 

As is well known, the gauge action Sg has Z3 symmetry, i.e., Sq is invariant when 
link variables on a time slice change as 



Ut ^ zUt 



where z is an element of Z3 group, 



z = e 



i<f> 



0, 



27r 47r 



(3-27) 
(3-28) 



3 ' 3 

But the fermion action is not invariant under the transformation. 

In the case of the imaginary chemical potential, the transformation of the 
link variables can be absorbed into ^ij in the fermion action. Conversely, the shift of 
the imaginary chemical potential. 



IJ^i ^ + 



(3-29) 
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Fig. 15. The equation of state correction Zi(p/r'*) vs. T/Tco for various /i/Tco by Bielefeld-Swansea 
collaboration. (Ref.EU)- 



can be managed as the change of the hnk variables, Ut zUt, which can be regarded 
as the simple change of the integration variables Ut to zUt- Therefore, the path 
integral gives the same result, 



' 1 ' 1 ' 1 ' 1 ' 1 ' 1 
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Fig. 16. Imaginary chemical potential study by de Forcrand and Philipsen. Critical Pc obtained 
with the imaginary chemical potential fxi is fitted by Pda/ii) = co + cl{afii)^ (left), and con- 
verted through Im/i = —iRe/i to obtain the critical lines in {T,jj,) (right). (Ref. 45,)). 



Under this constraint, we can perform the analytic continuation of the results 
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obtained with /x/ to the real //. The idea has been known for many years, and recently 
de Forcrand and Philipsen, and D'Elia and Lombardo have successfully performed 
simulations. In Ref . I45|) . the critical /3 is fitted as *^ 

/9 = oo + 02/^/ H (3-31) 

and transformed to the real chemical potential, 

/? = ao - aa/x^ + • • • (3-32) 

The result is shown in Fig |161 

Imaginary chemical potential may be considered as a special boundary condition. 

§4. Two Color 

Since SU(2) matrix U has the following property, 

U; = a2U^a2, (4-1) 

A{U,J^r = A{U\-fmU*)a2A{U,JrnU*)^2, (4-2) 

then {det A{x, x'; 7^)}* = det A{x, x'; 7*). 7* are a representation which satisfies the 
anti-commutation relations as 7^ and det A should not depend on the representation. 
**^ Consequently, 

det A{x, x'; 7^)* = det A{x, x'; 7^) (4-3) 

i.e., the fermion determinant is real. For reader's convenience, we show a compilation 
of color SU(2) finite density simulations in Table |31 

In Refs. IHTj) andl^S]). the authors have determined the lowest order, i.e., 0(i?^), 
effective Lagrangian including the chemical potential. 



where 



-FV^Tr {sb'^U'^B + BB^ - F^m^ReTr (^MS^ . (4-4) 



*' On the genuine phase transition hne, the partition function should be singular, but small 
finite lattice has only a relic of it. 

**' An explicit transformation between 7^ and 7^^ can be easily found: Since C'y^C — — 7jl, 
7m — V'y^V~^ with V = C75, where C — 27072 is the charge conjugation matrix and we use 
hermitian 7 matrices. 
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Fig. 17. The scaled chiral condensate, (tptl)} / {^ptp)o, (left) and number density, ub, (right) together 
with the chiral perturbation theory prediction. (Ref. EHJ) 



It has been shown in Ref. I69p that the numerical data are well described by the 
chiral perturbation theory. Fig llTI presents the chiral condensate and the number 
density, respectively, together with the leading order prediction, 

(V^V^) _ / 1 f \ _ / X < 1 

mo~\^ 'V8(^^)or'^" I 1(1-^) 

where 

X . ^ (4-8) 

Using extended QCD inequality, Kogut, Stephanov and Toublan have shown 
that in SU(2) color QCD at finite density, condensation must occur in scalar diquark 
channel.^''' The correlator of a meson, ■i/'F^/; is given by (TrG(2;, 0)/^G(0, x)r), where 
G = are quark propagators and taken as a matrix in Dirac and color space. At 
/i = 0, Ea. (|2-3p leads to = j^Aj^. From Cauchy-Schwarz inequality, one obtains 
Tr^^t < VTrAAWTvBBt Then 

TrG(x,0)rG(0,x)r = TrG(x, 0)r75G(x, 0)^75r 

< ^ Gix, 0)Gix, 0)t ^ r75G(x, 0)t75r (r75G(x, 0)t75r)^ 

= TiG{x,0)G{x,0)^ (4-9) 

The right-hand side corresponds to the pion propagator which behaves as G-j^ exp(— m 
at large x. Therefore 

C7e-"l^l < C7^e-"-l^l (4-10) 

This relationship holds at arbitrary large x only if m > m,r, i.e., the pion should be 
the lightest meson. This is the essence of QCD inequality. 

The above QCD inequality does not hold at finite density. For SU(2), however, 
there is another relationship, Va2Ga2V~^ , and using this the authors of Ref. I67|l 
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show that the mass of the scalar (O"*",/ = 0) diquark, ipc<^2lb''P^ is lowest, where 
ipc = ifFC. Therefore if diquark condensation occurs, it appears in this channel. 



Table II. Two Color 



Action 


parameters 


comments 


Ref. 


Plaq. + Wilson 


Nf =2, 1.4 ~ 2.0, 8^ X 2 


first calculation, pseudo-Fermion 
Method 


m 


Plaq. + KS 


P = 1.3, ruq = 0.05, 0.07, 6=* x 12 


symmetries and spectrum of lat- 
tice gauge theory 




Plaq. + KS 


(/3 = 1.3, rrig = 0.07,0.05), (/3 = 
1.5, m, = 0.1,0.07,0.05), (/3 = 
1.3 ~ 2.3, m, = 0.1), 6'' 


interquark potential 




Plaq. + adjoint KS 


Nf = 4, P ^ 2.0, nig = 0.01 - 0.1, 
4^ X 8 


adjoint dense matter, Two-Step 
Muhi-Boson algorithm 




Plaq. + KS 


Nf =8, 13 = 2.0, rUg = 0.025 ~ 0.2, 

4^ 6* 


fermion condensates 


60) 


Plaq. + KS 


Nf = 4, m, = 0.05, 8"*, 8^ X 4, 12=^ x 
6, 16* 


diquark condensation, tricritical 
point 




Plaq. + KS 


Nf = 4, /3 = 1.5, rrig = 0.1, 8'', 
12^ x 24 


phase diagram, diquark conden- 
sation 




Plaq. + Wilson 


Nf =2, 1.6, 4", 4^ X 8 


thermodynamical quantities, Linl 
by-Link update 


-64) 


Imp. + Wilson 


Nf =2,3, P^ 0.7 4-^, 4^ X 8,12, 


behavior of hadrons, Link-by- 
Link update 


65) 


Plaq. + KS fermions 


N/=8, 12 X 12 X 24 X 4, /? = 1.1, 
m,=0.05, 0.07, 0.1 


chiral condensation 


(Ml 



4.1. Diquark condensation - Super fluidity 



Wilson Line and Condensates on 12x12x12x6 lattice, chem pot=0.00, m=.05, ]amhda=0.005 



* ^ <qq> 



Wilson Line and Condetisates on 12x12x12x6 lattice, chem pot=0.30, m=.05, lainbda=0.Q05 
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Fig. 18. Wilson line and condensate vs. /3 calculated by Kogut, Toublan and Sinclair. Standard 
gauge -I- Staggered fermions. Nf = 4. Left : /j, = 0.0 Right : /j, = 0.3. (Ref. 163 



In Figs ll8[ we show the Wilson line, the chiral condensate (xx) the diquark 
condensate (xf2X) obtained by Kogut, Toublan and Sinclair Here x and x are 
staggered fermions. Using the staggered fermions, they simulated a system with a 



24 



S.Muroya, A.Nakamura, C.Nonaka and T.Takaishi 



diquark source,*^ 

S = xAx + jxcT2X^ + fx^^2X + Sg = ix, X^) ( _-r^T ) ( ) + 5g 

(4-11) 

After integrating out x and X; we obtain 

Z = y"pC/^detZ\Z\^ + |j|2e-^^G. (4-12) 

The fermionic measure is always positive. The source term acts as an external field 
in spin systems, and j should be extrapolated to zero at the end. Due to the square 
root, one cannot employ Hybrid Monte Carlo simulation, but the molecular-dynamics 
type algorithm is available to generate configurations. 

Fig ilHI (right) indicates the temperature behavior of the phase at fixed finite fj, 
with small j. When the temperature increases (the coupling increases), the Wilson 
line and the chiral condensate increase, which shows the confinement to deconfine- 
ment transition. The same behavior is seen at zero chemical potential shown in 
Fig |18l (left). A remarkable feature at finite n is that {x^ ct2x) has finite values at 
low temperature. Therefore at finite density and low temperature, there is a region 
in which the di-quark condensation occurs. From these observations, they suggest 
the phase shown in Fig |191 

The operator {x^(T2X) is & color singlet and cannot bring about color charge. 
Therefore this is not a color super conductivity, but it can be super fiuidity such 
as liquid ^He. This region belongs to the confinement phase, and this is probably 
because this diquark condensation phase is revealed by the color singlet external 
source. 

Schematic Phase Diagram of Diquark Condensation in the Temperature-Chemical Potential Plane 



Quark-Gluon Plasma 


<qq>=0 


2 


r <qq> HO 


<qq>=0 




Hadronic Matter 


Superfluid Phase 



Fig. 19. Schematic phase diagram in the (T, /i) plane by Kogut, Toublan and Sinclair. fRef. I63ll l 



Without the source terms, the expectation values of diquark states, (x^x)i ix^x) vanish. 
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4.2. Hadrons at finite density - In-medium effect ? 

In-medium hadron features are of great interest in high-energy phenomenology. 
Brown and Rho conjectured that 

^ = f . (4.13) 

which is now known as Brown Rho scahng, to explain the large low mass lepton 
pair enhancement observed in CERES Here m* and /* are the mass and the 
pion decay constant in the medium.'^ Based on the QCD sum rule, Hatsuda and 
Lee predicted a decrease of p meson mass as a function of //.^ Harada et al. 
showed that the vector meson mass vanishes at the critical density as a consequence 
of an effective theory with hidden local symmetry.^ Yokokawa et al. proposed 
simultaneous softening of a and p mesons associated with the chiral restoration.^ 
If lattice data show such a phenomenon, it is a new feature of QCD at finite density, 
and will become an interesting experimental target in future heavy ion experiments. 

Recently, Muroya, Nakamura and Nonaka reported a strange behavior of the 
propagator of a vector meson channel with chemical potential. They investigated 
finite density state of color SU(2) QCD with Wilson fermion. As for the gauge 
part of the lattice action, both plaquette and Iwasaki improved actions are adopted. 
In order to make a reliable update, the ratio of fermion determinants is calculated 
exactly through the Woodbury formula. As a result, the lattice is limited to a 
small size. They investigated the propagators of mesons and diquarks; pseudo-scalar 
(vr), scalar (oq), vector [p] mesons, and scalar diquark (qj^q), pseudo-scalar diquark 
(qlq) and axial-vector diquark {q^f_iq) are studied. At = 0, the propagators of the 
corresponding meson and diquark are degenerate with each other. With finite p, the 
propagators change their shapes. It is quite different from the finite temperature 
case, the propagator of the pseudo-scalar meson changes only slightly and the vector 
meson becomes lighter. The propagators of the corresponding diquarks show clear 
asymmetry in the chronological direction due to the charge corresponding to the 
chemical potential. However, the behavior of the obtained mass is very mild. Only 
the vector meson mass shows a marked change around pa = 0.7. This behavior 
might be the expected result of the phenomenological models. Figs l2UI are typical 
behaviors of masses of the pseudo-scalar, vector mesons, and di-quarks as a function 
of mu. 

The propagators of the scalar meson and pseudo scalar diquark decay so quickly 
that it is very difficult to catch signal against noise. 

Another possible explanation for the lightening of the mass of the mesonic chan- 
nel is a particle-hole pair around a Fermi surface. Hands, Kogut, Strouthos and Tran 
investigated the (2-|-l)-dimensional Gross-Neveu mode and discussed the effect of the 
Fermi surface. They evaluate the mesonic propagator in the high density state with a 
Fermi surface. In the mesonic channel, in addition to the ordinary quark-anti-quark 
pair which contributes as a massive meson, quark-hole gapless pair around Fermi 
surface may exist. In Ref.,^ authors investigated propagators with fixed momenta 
and found that at high density, a mesonic propagator whose momentum is up to 
some value shows massless behavior. On the other hand, propagators with larger 
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Fig. 20. Pseudo scalar and vector meson masses (left) and corresponding di-quark (baryon) masses 
(right). fRef.rflH. 

momenta behave massive. The existence of such a particle-hole pair is well known 
in the many-body theory and from direct results of the Fermi surface. In Ref.j^ 
a momentum dependent massless mode appears in the scalar channel and vector 
channel, but, in Ref. I74|l it appears in the vector channel only. In order to clarify 
the Fermi surface effect, we need larger lattice simulation and calculation of the a 
channel which also contains a disconnected diagram. 

§5. Related works 

5.1. Gross-Neveu model - Simple model for test 

QCD in four dimensions is highly nontrivial. A simple model may help us to 
understand the theory, and may give us a hint toward determining a correct route. 
The Gross-Neveu model has been studied for this purpose. 

Let us start from a model in two dimensions with Nambu-Jona-Lasinio type four 
fermion interaction which has a 'flavor' degree of freedom N. The action is given by 



TV 



5 = E 



k=l 



9 



where the Dirac matrices are 

71 = 72 = en, 75 = «7i72 = 0-3. 
When its mass m vanishes, Ea. ()5-1() possesses the chiral symmetry 



(5-1) 



(5-2) 



(5-3) 



We introduce two auxiliary fields, o"(x) and it{x), which correspond to —{g'^ /N)ipip 
and —{g'^/N)Tjjij^'tp, respectively, in the standard manner, 

^2 



9 
2N 
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N 
2? 



V^CTm^a* + m + a + i757r)V' + + ^^)- (5-4) 



Here we suppress the flavor index k. It is easy to verify that Ea. H5-4p is invariant 
under the following rotation together with Eq. l|5-|-{|) in a massless limit, 



cr \ ^ / cos 261 sin 261 \ ( cr 
vr / ^ I -sin26l cos26' / I vr 



(5-5) 



In a large limit, the fields a and vr can be replaced by their mean fields, 
and we may rotate them so that tt = 0. Then we have a model with only the first 
four-fermion interaction term in Eq. (|5-H) . or 



S = ij{-i^d^ + a)ij + ^a\ (5-6) 



Although the model has no confinement, it has the following features common 
to QCD: 

• it is an asymptotically free theory, 

• it shows spontaneous breakdown of a (discrete) chiral symmetry, i.e., ij) — > 75^. 
The Gross-Neveu model is well studied in the continuum, 

• the phase structure at finite temperature and density is obtained in A'^ — > 00 
limit by a mean field approach. 

• Finite N case can be treated by the factorized S-matrix method.^ 

• There is a conjecture that a generalization of the Gross-Neveu model from U{N) 
to 0{N) flavor symmetry leads to the appearance of a pairing condensate at 
high density.!^ 

There are two lattice studies of the Gross-Neveu model in the literature. Karsch, 
Kogut and Wyld analysed the staggered fermion version of Ea. H5-4p .l^ 

^ 1 
S = Y.x'^x' + ^a^ (5-7) 

k=i ^ 
where the fermion matrix is given by 

(5-8) 

All elements of the fermion matrix are real, therefore its determinant is also real. 
In Ref. I82|) . the Langevin algorithm was used to simulate the system with N = 12. 
Figl^ shows typical simulation results. 

Izubuchi, Noaki and Ukawa investigated a lattice version of Eg 15- II using Wilson 
fermions at finite temperature and density,^Si 

S = - ■^['il^{x){r - "fi)i'{x + i)V5(x)(r -|- 'Ji)'il^{x - 1)] 
la 
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Fig. 21. Left : (a) as a function of ^ on latticetj |of size 60 x Nr,Nr = (\,8, 1 0, 14, 20. Right 
Coexisting states at Hc- Ref i82i . 



— [^{x){r - 72)e~^'*V(2; + 2)^{x){r + 7i)e^Xx - 1)] 



+ -{2r + 5m)i}{x)i;{x) - ^{i;{x)ijix)f - ^{i;{x)i^,^{x)f . (5-9) 



Here they introduce two couplings and to controle the continuum hmit. The 
lattice effective potential is constructed with two auxiliary fields, a and vr as in the 
continuum, 



Vr 



y g2 + 2Ci) + - 2Co) vri + f-^ - Co + 2C2 ] al 



9l 



+i-(ai + vri)log^l±^ + 0(a3), (5-10) 



where Co, Ci and C2 are constant. In order to get the correct chiral symmetric 
theory in the continuum, we tune three parameters, g-,^, g^ and 5mi so that 



\ = \ + AC2 + 0{a) 
9i 9i 



6m L 



-2Ci + 0{a^). 



(5-11) 



5.2. NJL model 



At present, there are no reliable methods for evaluating low temperature and 
high density regions in nonperturbative QCD calculation. Therefore effective the- 
ories such as Nambu-Jona-Lashino (NJL) model play an important role and the 
comparison of the results of lattice QCD for a simple case to the results of effec- 
tive theories is a meaningful approach in order to understand the superconducting 
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behavior at high density. Hands and Walters apply the lattice analyses to the NJL 
model in (3+l)-dimension, focusing on BCS diquark condensation-SS) *) The action 
of this model is written by 

5 = V'f-^V' + 4y"(f^^ + vf -vf), (5-12) 

^ X 

where the bispinor is written in terms of staggered isospinor fermions fields, ip^"^ = 
(X)X*'^) s-iid a and the triplet vf are real auxiliary fields defined on the dual site x. 
In the fermion matrix A they introduce the diquark source term j, 

where the matrix M is defined by 

+ ) + 271206 xy ] 

l/=l,2 

+—Sxy V {a{x)6P'^ + ie{x)Tf{x) ■ r^^). (5-14) 
16 

Here mo is the bare mass and r]u{x), e(x) are the phases (|_i)^'o+ --a;i/-i g^j^^j ^_-^-^xo+xi+'. 
respectively. The r are Pauli matrices and {x, x) represents the set of 16 dual lattice 
sites neighboring x. They measure diquark condensate {qq) as a function of chemical 



0.7 - (xx) I— e— I 



riB l—B- 




Fig. 22. The diquark condensate, (qq) as a function 
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J 

of /i(left). (qq) as a function of j (right). 



potential, extrapolating to j 0. From this analysis, they determine the existence 
of a nonzero BCS condensate at high baryon density (Fig. |221(left)) . However, the 
behavior of {qq) as a function of fi depends on the evaluation of the data at j 0. 
They ignore the values of {qq) at small j, because at high /i, {qq) at small j falls 
rapidly away from the tendency of {qq) at large j (Fig. EHright)). They discuss that 

*•* The (2+l)-dimensional lattice NJL model at finite density has been simulated, but BCS 
condensation is not observed!^ 
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the low-j discrepancy comes from the finite volume effect and the partially quenched 
approximation and some issues concerning thermodynamic limit need to be resolved 
for a definite conclusion. 

5.3. Strong coupling calculation 

It is comfortable if we have any analytical or semi-analytical method to study 
QCD at finite baryon density in non-perturbative region. Strong coupling approxi- 
mation of lattice QCD is one of such possibilities. Indeed the lattice formulation was 
introduced to allow the strong coupling expansion and the confinement was shown 
there. We can also compare our numerical data to the analytical calculations and 
take a deep view of lattice results. 

When the gauge coupling g is very large, pSc = ^/g^Sc may drop, 



VUVi/jViP e-^""^ . (5-15) 

In this case, the integrations over time like link variables, U4{x), and space-like ones, 
Ui{x), are decouples, and one can integrate over the gauge fields. 

The first strong coupling calculation was performed by Ilgenfritz and Kripf- 
gan^^ and Damgaard, Hochberg and Kawamoto, together with the mean field 
approximation. Bilic, Demeterfi and Petersson extended the calculation to the next 
order of 1/g'^W^ Recently Nishida, Fukushima and Hatsuda applied the strong cou- 
pling method to two-color QCD to study the dynamics of the theory.!^ 

For warming-up, let us first study a very simple case in a pedagogical way: 

Z= [ dUdipid'iljidip2d'ilJ2e^'^^'-'^'^^'^\ (5-16) 



i.e., two fermion fields, ipf, connected by a link filed U, where a = 1, 2, • • • ,Nc 

are color index. Series expansion of the exponential is finie due to the nil potency of 
the Grassmann variable. With the help of group integration properties. 



/ 



dU 



1 



dUUa,bU^c,d = :^Sa,dSb,c = d 



a ' b 



- c 



a, ► b, 

1 az ► bz 

dUUa^^biUa2,b2Ua3,b3 = ;^eaia2a3^feif'263 = ^ 
J" dUU^ ai,biU^ a2,b2^^ a3,bi, — ^^aict2Ci3^feib2''3 (^'^''') 



where Ua^ is (a, b) element of the matrix U, we obtain 

Z = l + \miM2 + {-B1B2 + B2B1) + 



^^MiM2+i-BiB2+B2Bi)^ j^g.-j^g^ 
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where 



B = —f i/;"ii/;°2,/,.a3 

— 21 '^010203 Y^i r« Yl 



(5-19) 



More systematic and precise formulae are given in Ref.® ' Using the identities with 
auxihary bosonic field (j) and fermionic field 6, 



(5-20) 
(5-21) 



we linearize M and B fields, where M =* (M1M2). 

Now we come back to the realistic case. We integrate first the spatial link 
variables for the staggered fermions, 



X{x) 



-fir 



2r 



x,y+4 



x,y 



(y) - x{x)x{x)x{y)x{y) 



{x,y) 



32N^{N, - 1) E[x(^)x(^)]'[x(^)x(y)]^ i (5-22) 



where an anisotropy parameter of time-like and space- like lattice spacing, r = at/ag, 
is introduced to controle the temperature. Introducing auxiliary fields. A, we can 
linearize the exponential. 



VUiVxDxVX exp 



x) — m)V ^{x,y){X{y) — m) 



x,y 



X{x) 



1 



-fir 



x,y+A 



xiy) } (5-23) 



where d is the spatial dimension (= 3), and the matrix V{x,y) is defined by 



1 



(5-24) 



k=l 



The integration over x and x can be performed. Then we can further integrate 
over C/4. When A is constant (mean field approximation), then 



Z = 2 cosh{NtNcfir) + 



sinh(iVc + l)NtX' 
sinh(A^tA') 



(5-25) 
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where A' = sinh^^[(A + m)r]. A is determined by minimizing the free energy F = 
NtNc^? /d — logZ, i.e., dF/dX = 0. The chiral condensation is given by the A by 
ixx) = 2NcX/d. In the infinite couphng hmit, the transition value of the chemical 
potential is given as 

//o = -sinh-nAor)-^, (5-26) 
r dr 

where 

Ao = -i-(Vl + d2r4 - 1)1/2. (5.27) 
rv2 

The transition point with l/g^ correction can be found in 




Fig. 23. Strong coupling calculations by Bilic, Demeterfi and Petersson (Ref. I88II '). Chiral con- 
densate as a function of /x for various temperatures at &/g^ = 2. The dashed line represents 
the instability region, (left) Comparison of the strong coupling calculation of the quark number 
density with numerical data for A'^t = 4 and &/g^ = 4.8. The points are data of Re£^''' (right) 

In Fig |23l we show the chiral condensate as a function of the chemical potential 
for various temperatures at 6/5^ = 2 taken from Ref.l^ where the 1/5^ correction 
is included. They calculated also the quark number density as a function of the 
chemical potential and found good agreement with the numerical calculation given 
by Hasenfratz and Toussaint.'^^ 

Karsch and Mutter have given an elegant form, 'monomer-dimer-polymer' sys- 
tem.l^ For the fermion matrix of the form, 

Z\ = 2ma + ^C/(x,y)C(x,y), (5-28) 

{xy) 

Z = J PVi^?V'e2'""S^^(")i7(.,)F(x,y). (5-29) 

where 
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= 1 - ^p{x,y)Mix)Miy) + y)M(x)M(y))2 - ^{pix,y)M{x)M{y)f 

+C(x, yfB{x)B{y) + ({y, xfB{y)B{x). (5-30) 

Here M{x), B{x) and B{x) are 'meson', 'baryon' and 'antibaryon' fields. For stag- 
gered fermions, the quark field tp{x) has only color degrees of freedom, not Dirac 
indices. Meson, baryon and antibaryon fields are given as, 

M{x) = i}{x)'^{x) = ^ = ^ Ma{x), (5-31) 

a=l,2,3 a=l,2,3 

B{x) = ip^{x)ip'^{x)4^^{x), 

B{x) = i;^{x)'ii;'^{x)i;^{x), (5-32) 

respectively, and 

^2maM{x) ^ ^ ^ 2ma {Mi{x) + M2(x) + Msix)) 

+ {2maf {Mi{x)M2{x) + M2{x)M3{x) + Ms{x)Mi{x)) 

+ {2maf {Mi{x)M2{x)M3{x)) . (5-33) 

The last equation gives 'monomer' contributions. p{x, y) is the 'mesonic' link factor, 

p{x, y) = C,{x, y)C{y, x), (5-34) 

from which three types of "dimers", p{x,y), p{x,y)'^ and p{x,y)^, are constructed. 
For staggered fermions, 

±i: • f::4tM^= 1.2.3) <^-^^) 

with r] being the staggered sign factor. 

Integration over Grassmann variables, ijj and ip, is straightforward. Nonvanishing 
contributions appear only if each site of the lattice is occupied either by three mesons, 
Afi(;r)Af2(x)M3(x) or by a baryon-antibaryon pair, B{x)B{x). Finally, the partition 
function is written as a sum over monomer-dimer loop configurations K, 



Z{2ma, as/at, p.at) = ^ wk, (5-36) 



K 

where 

2(# of Dimer lines in t direction) 



WK = {2map of Monomer) (^^J 



-j^ ^ (# of the first and second type Dimers) 
3 

w{x) and w{C) are site and loop weights associated with the configuration 



X ( - 1 n^w{x)ncw{C). (5-37) 
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Fig. 24. (a) The chir^l condensate (ipip) — V^[\6iog Z / d{2'ma) and baryon number density na^ 
{3VNt)^^dlogZ/d{fiat). (b) The average sign of the Boltzmann weights. (Ref. 91)) 



In Ea. ()5-37|) . wk is not always positive. A sign ax for each monomer-dimer loop 
configuration is introduced as 

Wk = ctkIwkI, (5-38) 
and the simulation was carried out with the positive Boltzmann weight \wk\, 

(O) = (5-39) 

Fig |24l shows the chiral condensate and baryon number density as a function of fia 
for ma = 0.1 on 8^ x 4 lattice obtained by Karsch and Mutter. They determined 
that /icO = 0.69 lb 0.015. The average of the sign factor drops at the critical /j. shown 
in FigEKb). 

5.4. Canonical ensemble 

A thermodynamical system can be studied either as a grand canonical formu- 
lation or as a canonical one. When the grand canonical partition function Z{iJ,) 
is expanded in terms of the fugacity, = exp(/i/T), each coefficient is a canonical 
partition function Zjy with a fixed quark baryon number N, 

Z{f,/T)= ^ (e^) Z^. (5-40) 

^=-00 
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The canonical partition function can be formally obtained from the grand canonical 
partition function with imaginary chemical 

potentiaipm 



2-K 



r d(t>Z{i4> ^ fi/T) = V C(n, n)C"C", (5-41) 

•^0 n,n 



where C = e*''',^* = e"**^. The canonical partition function with the quark baryon 
number N = n — n has contributions from (^"(^*" terms in the grand partition func- 
tion. Miller and Redlich developed a general lattice canonical formulation, and 
investigated it using the hopping parameter expansion. 

To date, only one numerical simulation was performed by Engels et al.!^ They 
used Wilson fermions, Ea. ()2-7() . and the expression (|2-17() . Then the coefficient of 
^"C*" should include Therefore for a heavy quark system, i.e., for small k, 

the leading contribution comes from the n = sector. An explicit expression of the 
expansion is very complicated, and yet they succeeded in calculating up to = 12 
in the quenched approximation. The baryon number density is given as 



S/3 



(5-42) 



iV = 12 on 8^ X 2 lattice at T ~ Tc = 270 MeV corresponds to ub - 0.15 fm^^ , i.e., 
approximately nuclear matter density. 

The heavy quark potential shows finite values for the large distance shown in 
Fig |251 = 6 on this lattice at T ~ 0.86 Tc corresponds to = 0.049 fm~^. The 
heavy quark potential is screened even below the nuclear matter density. 
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Fig. 25. Heavy quark potential for T ~ QMTc and B = and 6. on 16^ x 4. 



5.5. Density of states method 

In lattice QCD, the average value of an operator (O) is given by 



1 



VUO[U] dQiA{U)e 



(5-43) 
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with 

Z = I VU det A{U)e~f^^^^^\ (5-44) 



It is impractical to directly perform this integration of very high dimensions. Instead, 
we can use a importance sampling where configm'ations Ui are generated with the 
probability distribution: 

P{U) oc det Z\([/)e-^^3(^) (5-45) 
1 ^ 

and (O) is given by lira — 0[Ui] where N is the number of configurations. 

i 

If one knows the density of states of some parameters, then the integration can 
be reduced to an integral of few dimensions. Let us assume a one-parameter case 
and define the density of states as 

p{E) = [ dU6{E -x{U))g{U) (5-46) 



where £" is a certain parameter associated with x{U) for which we construct the 
density of states p{E), and g{U) is any function of U. Then using the density of 
states p{E), (O) is given by 

^0) = Y j dEp{E){OdeiAe~''^^/g(U))E, (5-47) 

Z^ = I dEp{E) {det Ae-^^^/g{U))E, (5-48) 
where (O) e stands for the micro canonical average with fixed E, defined as 

{0)e = ^ / dU6iE - x{U))g{U)0{U). (5-49) 

In this way one can define any p{E) which depends on x{U) and g{U)^ and (O) 
is rewritten by using the density of states p{E). In the following, we show three 
examples of x{U) and g{U) which have been actually studied in the literature. 

5.5.1. Density of states of energy 

The simplest choice of the parameter for the density of states may be the energy 
of the gauge action. In this case we set 

x{U) = Sg{U), (5-50) 

g{U) = 1. (5-51) 
Substituting these equations to Eas. ()5-46() - ()5-49|) we obtain 

p{E) = j dU5{E - Sg), (5-52) 
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{0) = ^ [ dEp{E)e-^^{OdetA)E, (5-53) 

Zp = j dEp{E)e-'^^ {det A) E. (5-54) 

This method including dynamical fermion was applied in compact QEEC^ and 
then in QCdP Both were simulated at zero density. Considering the complex phase 
in the above equations, one may apply this method for the finite density case. If we 
explicitly introduce the complex phase in the equations, Eas. (|5-53|) and (|5-54p are 
rewritten as 

(0) = ^ / dEp{E)e-'^^{0\detA\e''^)E, (5-55) 
Zp J 

Zp = j dEp{E)e-'^^{\detA\e''^)E. (5-56) 

Since there is still the sign problem, e.g. e*^ fluctuates, it is not obvious that this 
works better than the standard simulation. If we omit the phase factor, resulting in 
simulating the finite isospin model, the method should work. Fig |26l shows a com- 
parison of < 'i/'V' > between results from the density of states method and those from 
the standard Monte Carlo simulation(e.g. importance sampling) The density of 
states method reproduce the results from the standard Monte Carlo simulation well. 



1.5 r 
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□ Density of States 
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Fig. 26. < tp''P > for ^^/ = 4 on a 4** lattice at /j / = 0.2 and niq — 0.025. The standard Monte 
Carlo results were obtained by the R-algorithmP^ Ref ll()5|l 
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5.5.2. Density of states of phase 

GocksclPJ used the complex phase of the determinant as the parameter for the 
density of states and studied the chiral condensate {^^) and the phase diagram of 
QCD at infinite gauge couphng. In Ref. IM)) he sets 

x{U)=9{U), (5-57) 
g{U) = \detA{fi)\e-'^'^o. (5-58) 



Then we obtain 



p{E)= j dU5{E -e{U))\detA{n)\e-f^^\ (5-59) 
(O) = ^ y" dEp{E)e'''{0)E, (5-60) 



'' JE 



Zp = j dEp{E)e"^. (5-61) 

In order to obtain p(-E'), simulations were carried out according to the probability 
P{U)dU ~ \ deiA{p)\e-^^^dU and the histo gram of E, which forms the density 
of state p{E) in the end, was constructed. In the actual simulations, in order to 
speed up the process, the period of E was divided into bins and the simulations were 
performed for each bin. Each bin has overlaps with the neighboring bins to match 
the adjacent density of states. The updating process is constrained in the range 
of each bin. Therefore the configuration leading out of the range is not accepted. 
Nevertheless the event has to be recorded in order to form the correct histogram .1^^ 

5.5.3. Density of states of observable 

One may use observables to form the density of states. In Refs. I1U1|) ancP^ 
the random matrix theory for finite density QCD was studied by making use of the 
density of states for the quark number density. Originally they call the method 
"factorization method" which is considered as a variety of the density of states 
method. 

Setting 

x{U) = 0{U), (5-62) 

and 

g{U) = \deiA{p)\e-^^\ (5-63) 



we obtain 



p{E)= J dU6{E -0{U))\detA{p)\e-^-'^\ (5-64) 
(0) = yJ dEp{E)E{e'')E, (5-65) 



Zp= I dEp{E){e'^)E. (5-66) 
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If the operator O takes a complex number, we may evaluate real and complex 
terms independently. First we rewrite (O) as 

(O) = {Or) + i{Oi) (5-67) 

where Or and Oj stand for Re(0) and lm(0) respectively. Then, defining the density 
of states Pr{E) and pi{E) as 

Pr{E) = J dU6iE-OR)\detA{i^)\e-^^\ (5-68) 

pi{E) = j dU6{E - Oi)\detA{p)\e-'^^\ (5-69) 
we obtain (Or) and (Oj) as given bellow. For i = R, I, 

(Oi) = dEpi,{E)E{e'')i,E. (5-70) 

Zi = J dEp.,{E){e%,E, (5-71) 

where 

{e%,E = [ dU6{E - 0,)e''\ det A{p)\e-'''^ . (5-72) 
Finally the real part of (O) is given by 

Re(0) = {Or) - {6i) (5-73) 

where 

{Or) = dEpR{E)E{cose)R^E (5-74) 

{Oi) = ^ J dEpi{E)E{sm 9)j,e (5-75) 

and C = Zr = Zj. In Ref.,!^^ to obtain {cos 9) r^e and (sin6')/^£;, they performed 
simulations with the constrained partition functions which are actually the density 
of states as given in Ea. (|5-68|) and (|5-69j) . The density of states itself was also 
evaluated through those simulations. They argue that since the simulations are 
performed with the constrained partition function which forces the simulations to 
sample the important region, the results are obtained accurately. 

5.6. Complex Langevin 

Stochastic quantization by Parisi and Wu does not rely on the path integral, and 
its simulation is realized by the Langevin equation 

Pl ParisPZJ and 

KlaudePHJ 

proposed independently an extension of the Langevin algorithm to a system with a 
complex action. 
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For the dynamical variable x which is subject to the (Euclidean) action S{x), 
let us consider a stochastic process along the "fictitious" time r, 

dxlr) dS , , ,^ 
-P = -— + r?(r), (5-76) 
OT dx 

where r;(r) is a Gaussian white noise characterized by {rj^r)) = and {viT)v{'^')) — 
26 (t — t'). From the Langevin equation, we can derive an equivalent Fokker-Planck 
equation, 

d , d / d dS\ 



By virtue of the positive semi-definite property of the Fokker-Planck Hamiltonian, 
Hpp = — + we can obtain a stationary solution as, 

P(X, t) >Peq(x) OC C^^^^^. (5-78) 

t^oo 

Therefore, putting the solution of the Langevin equation into the quantities 0{x) 
and taking the average over the ensemble of the random noise r], we can obtain the 
expectation value (0(x)) with probability measure Peg (■^) ^ ^ ■ This is a rough 
sketch of the basic procedure of the stochastic quantization. In practical use, instead 
of an ensemble average, the long time average is taken over the simulation time r 
after thermalization, where we assume ergodicity. 

Once the action S{x) becomes complex, the situation is not trivial. Putting 
a complex action into the above Langevin equation leads to the complex-valued 
extended version, 

^ = -^^ + ,(r) (5-79) 
ar dz 

with a complex variable z{t) = x(r) -|- iy{T) and real noise ?/(t). The corresponding 
Fokker-Planck equation is defined for a complex distribution function Pc{x) of the 
real variable x which satisfies, 

{0{z{T)))rj = J dxPc{x)0{x). (5-80) 

In addition to the problem of the interpretation of the complex generalized 
distribution function Pc{x), in this case, Fokker-Planck Hamiltonian Hpp loses the 
positive semi-definite property and there exists no general proof for the convergence. 

However, if the spectrum of the Fokker-Planck Hamiltonian is composed of very 
small imaginary part and positive semi-definite real part, it is expected that a com- 
plete set Zn with positive semi-definite eigenvalue En and stationary solution e"'^^^-' 
exist. Fo r a simple case, this conjecture has been checked explicitly by Klauder and 
Petersen 

Karsch and Wyld applied the algorithm to a three dimensional SU(3) spin model 
with the chemical potential.^^ The Langevin simulations converge even for large 
values of the chemical potential, and are in agreement with exact solutions at the 
strong coupling limit. 
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Unfortunately the complex Langevin method sometimes gives incorrect results. 
For some simple noncompact variable case, the Parisi and Klauder algorithm can 
be proved to work well if the partition function satisfies certain conditions of the 
analytic behavior, but for the compact variable case, it is difficult to judge if the 
algorithm gives correct results or not. 

Extensive use of the fictitious time as a redundant variable introducing a ker- 
nel in the Langevin equation has been discussed by the Waseda and Siegen group 
intensively.^^ A kernel K{z) is introduced as, 

dzir) , , dS dK , , , ^ 

where ^(r) is a noise satisfying (^(r)) = and (^(T)^(r')) = 2K{z)5{t — r'). The 
corresponding Fokker-Planck equation becomes, 

^Pc{x; t) = -HfpPc{x; t), (5-82) 

where 

d , , f d dS\ , , 

Therefore, the role of the kernel is to change relaxation process only keeping the 
same asymptotic behavior, 

Pc{x,t) >Peg{x) oc e-^("). (5-84) 

Thereby, a choice of the appropriate kernel may provide a better relaxation pro- 
cess without changing the physical result. Indeed in some cases, the use of a ker- 
nel improves the relaxation process and converges the random walk. For example, 
Okamoto et al. applied a constant kernel method to the polynomial model and 
succeeded in ext end ing Klauder and Petersen's convergent region in the complex 
eigenvalue space.^^ Though a pole in the second Rieman sheet sometimes leads 
the constant kernel method to an incorrect result, an appropriate choice of the field 
dependent kernel solves the problem in the case of the polynomial model.^^ For the 
simple models, the origin of the problems in the complex Langevin method is well 
analyzed and some answers with kerneled Langevin equation exist (see more details 
in Ref. 1115]) ^. However, the general recipes for the complex system are still unclear 
and how to manage the finite density QCD is still open. 

5.7. Finite isospin model 

Imagine a system in which the chemical potentials of u and d quarks have the 
opposite signs to each other, 

fJ-u = -fJ-d = M (5-85) 

The system is called a finite isospin model^^^- or the isovector model.!^ 
The fermion determinant of this system is 

det A{fiu) det A{fj,d) = det Z\(/i) det Z\(-/u) = | det Z\(/u)|^ (5-86) 
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where we use the relationship (|2-3|) . Therefore the finite isospin model is considered 
as a two-flavor system where one drops the phase of the fermion determinant, and 
of course Monte Carlo calculation is possible. 

Son and Stephanov have obtained the phase of the system shown in Figl^Tp^ 
and the qualitative features of the diagram were confirmed by numerical simulation 
by Kogut and Sinclair.EHl See Figs. l28l 
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Fig. 27. Finite Iso-spin phase diagram calculated by Son and Stephanov. fRef. IllSl ). 
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Fig. 28. Numerical results of finite iso-spin model by Kogut and Sinclair. Left : i{psi'y^T2^) at 
/i/ = 0.8 has finite value when the temperature is low, and decreases as T increases. Right : 
Wilson line at /i/ = 0.8 at low T is zero and increases when T increases. (Ref. [TT7|). 



§6. Summary 



After the brief introduction of the lattice QCD actions with the chemical poten- 
tial and their characters, we have seen many activities in the field. Recent progress 



Lattice QCD at Finite Density 



43 



in small /i and finite T regions is very promising. Indeed until several years ago, we 
did not expect that any physically relevant data were provided by lattice QCD. This 
is the area which recent high energy heavy ion experiments are investigating. We 
anticipate now lattice data in the next several years which can be compared with 
experimental results. For that algorithm development is important. 

Still it is difficult to study finite density QCD in large baryon density and low 
temperature regions, where many interesting phenomena such as color super condic- 
tivity and partial chiral restoration may occur. Two color QCD provides us hints 
on what may occur. Lattice simulations discussed in Sec|l]show us many interest- 
ing phenomena. Some of them are relevant for QCD and some not. In any case, 
they are interesting features in field theories and it is desirable to investigate their 
mechanisms. 

We discussed many attempts in SecEl related with lattice field theory at finite 
density, i.e., 2-dimensional model, lattice NJL model, application of the strong cou- 
pling technique, density of state method, complex Langevin and finite isospin model. 
Surely there are more in the literature. We have made a survey of these works by 
hoping that a new idea appears from these studies, i.e., taking a lesson from the 
past. 
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Appendix A 

Gibbs formula to condense the fermion determinant 

In lattice QCD simulations of finite baryon density system, we often need to 
calculate the fermion determinant. Gibbs has obtained a useful expression for this 
purpose. Here we describe a brief derivation of the formula following 55), 37). Our 
proof given below is not very smart, but we hope that it helps the readers who need 
the formula. 

We consider KS fermions. H2-8() . and employ Eas. H2- 16(1 . (|2- 17(1 . *^ The fermion 
matrix Z\ is (3 x NxNyN^Nt)'^ complex matrix. We take the axial gauge as Ut{x) = 1 



We renormalize a factor two in 12-8II . 
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except at t = Nf — 1, the fermion matrix can be written in time slices as follows: 



A 



Bo 


1 


••• 







-1 


Bi 


1 











-1 


B2 ••• 


1 











... _i 


BNt-2 


1 










-1 


Bwt-i 



(A-1) 



where each column and row corresponds to the time slice t = 0,1,2, - ■■ Nt — 1 

and Ujy^^i is the abbriviation of a {NcVg) x {NcVg) diagonal matrix whose clemets 
are Ut{x,Nt — l)6g^gi. Here we use the anti-periodic boundary condition along t 
direction. The element of the above matrix is a (NcVg) x [NcVg) block matrix, where 
Vs = N^NyNz- Bi stands for the spatial derivative and quark mass parts of the 
matrix at i-th time slice (i = 0, • • • , A^^ — 1). The chemical potential dependence has 
been summarized to the last time slice only. Let us denote U^^^ le^*/^ as ^ for the 
abbreviation. 

Let us define A by multiplying S only to the last column of A, 

\ 



( Bo 1 
-1 Bi 
-1 



A 



1 





1 

BNt-2 
-1 



B 



V -H 
We assume that Nt is even. Since 

det(H) = det(C/jv,-ie^*'^) = e^^'^^^". 



(A-2) 



(A-3) 



the determinant of A is related to detZ\ as detZ\ 
following matrix whose determinant is one, 

/ 1 5o 
10 

1 B2 

1 

B = 



V 



det Ae^^'^*'^ . Multiplying the 





BNt- 

1 



to the fermion matrix A from the left, we obtain. 



BxA 
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1 + BqBi Bo 

Bi 1 

-1 I + B2B3 B2 

-1 S3 1 



-BNt-2^ 



V 



1 





(1 + i?Ar^_2-BjVt-i)S' 
— 1 BiVt-iH 



Permuting the first column to the last and multiplying -1 to the last column, det A 
is rewritten as, 



det 



/ 1 + BoBi 


Bo 





1 


\ 


Bi 


1 








1 


-1 





1 + B2B3 B2 










-1 


B3 1 












(1 + Bf^^_2BNt-i)^ 


BNt-2^ 


V 






-1 Bn^-i^ 


^ J 



Denoting as, 
and 



BjBj+i + 1 Bj 



B 



Bj 1 
1 



Bj+i 1 
1 



I2 



1 
1 



(A-4) 
(A-5) 



the determinant of the fermion matrix can be written as, 

/ f^oi I2 

-I2 ^23 

— 12 i745 



det(Z\) = det 



V 



(A-6) 



-I2 f2Nt-2,Nt-l^ J 



Note that QNt-2,Nt-i^ is (3^^ x 2) times {'iVg x 2) matrix. 

We multiply a matrix of which determinant is 1 from the right. 



— I2 J?23 



det (4) = det 



— 12 i?45 



/ a 



det 



01 

-I9 Q- 



\ 



-I2 QNt-2,Nt-l'^ I 





/ 12 



12 



12 / 



23 



— 12 i?45 



-I2 QNt-2,Nt-\^ I 
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\ 



det(i7oi) det 



— 12 i?45 

V 



-I2 f2Nt-2,Nt-l^ I 



det (J7oi^23 • • • ■^^Aft-e.Aft-s) 

— 1 /7Arj_2,Art-l" + (^01 ' ' ' ^Nt-Z^Nt-"!)'^ 



where we used a formula for bfocked matrixes Aif 



Axx 





A2n 



A 



22 



A 



2n 



A. 



n2 



An 



Taking the similar step recursively, we can reduced the matrix up to two-time 
slices. We obtain finally, 



det(Z\) 



detzi X e-'i'^^^ti^ 
det(l + e^*^P) X 

det {e^ti^(e-^ti^ + i^)} x e-^^-^'/^ 

det [P + e-^*'^) X (e3^-^*A')2 x g-^^-^*/^ 

det (P + e"^*/^) X e3^=^*^. 



(A-9) 



where 



P 



n 



Nt-2 



i=0,2,4''^ij+l 



Nt-1 



(A-10) 



The rank of the matrix P + e~^*^ is 3 x {2Vs). Therefore if Nt > 2, we obtain 
a condensed smaller matrix. To estimate the determinant, LU decomposition is 
fastest for small size lattices on many computers if there is enough memory. If 
we have eigenvalues of the matrix P, {Aj}, which do not depend on the chemical 
potential /i, we obtain Ea. ()3-2ip . which gives us values of the fermion determinant 
at any ^. 

Appendix B 

Meson and haryon propagators in color SU(2) space 



A q — q bound state belonging to 1 representation in 2 2 is called a meson. 
With a quark field ip{x), the meson operators are constructed as 



(A-7) 
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with r being Dirac's 7 matrices, i"' = 1 corresponds to a scalar meson, -T = 75 
corresponds to a pseudoscalar meson, i"' = 7^ corresponds to a vector meson, and 
r = 'j^'jfj, to an axial vector meson. 

In SU(2), 1 representation also exists in 2 2 and q — q bound state can be 
realized as a color singlet state. It means that a baryon in SU(2) theory is a color 
singlet di-quark state. On this point, color-SU(2) and color-SU(3) are essentially 
different. The di-quark operator is given as, 

with a and b being color indices, f stands for Pauli matrix T2 or tt2 in the flavor 
SU(2) space. In order to satisfy the fermionic property of the quark, the di-quark 
state must be totally antisymmetric utilizing flavor degrees of freedom. We must 
take iso-singlet channel, f = T2, for i"' = 1, 75 and 7/^75 • Iso- vector di-quarks are 
realized for the vector channel, F = 7^. Because of the charge conjugation C, we 
must turn our attention to the parity of the di-quark states; F = 1 corresponds to a 
pseudoscalar di-quark, F = 75 corresponds to a scalar di-quark, F = corresponds 
to an axial vector di-quark, and F = 757^ to a vector di-quark. 

By virtue of the simple structure of SU(2), propagators of meson and di-quark 
are degenerate in a vacuum state, i.e. , propagators of pseudoscalar meson and scalar 
di-quark, and vector meson and axial vector di-quark are degenerate. 
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Fig. 29. Hadron propagators at /xa = (left) and jj-a — 0.8 (right). Pseudoscalar and vector 
mesons and scalar and axial vector baryons are shown. 



Figure H29|) shows the propagators on a 4^ x 8 lattice with spatially periodic 
boundary condition at /3 = 0.7 and k = 0.150. At = (|29k ). pseudoscalar meson 
and scalar di-quark, and vector boson and axial vector di-quark are degenerate. With 
the finite the fiB term in the W makes the di-quark propagators asymmetric in the 
time-direction. On the other hand, the meson propagators keep the symmetric shape. 
Therefore, degeneracy between the meson propagator and di-quark propagators is 
solved with finite 

The asymmetry of the di-quark propagator is caused by the difference of the 
charge of " di-quark" and " anti-di-quark" . The chronological propagation is provided 
by a particle with charge Q and of which the chemical potential is /iQ. On the other 
hand, the anti-chronological propagation is governed by an anti-particle with charge 
—Q and which is affected by the chemical potential as —fiQ. Taking this asymmetry 
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into account, we adopt 

as a fitting function to the di-quark propagator, where Ci and C2 must keep the 
constraint on the bosonic boundary condition of di-quark. For the mesons, even in 
finite density we can use a usual form as, 

G^t = C7ie-™* + C72e-™(^^-*). 

Appendix C 

Thermodynamical variables on the lattice 



The path integral form of the partition function is given aJ^ 

Z = Tr(e-^(^-^^)) 



g-/3(£-/.n)y 



(C-1) 



The energy density e, and the number density n are given by 

^ f d /i 9 \ ^, 



On the lattice, (5 = Ntat, and 

d _ I d 

dp ~ Wtd^t 



(C-3) 



We must be careful that the temperature, T, changes when we vary the lattice 
spacing at which is a function of the coupling constant. 

A simple and safe way is to write down a formula in terms of non-dimensional 
quantities and to start to calculate. For example, the number density is given in the 
continuum world as 



We may write a formula 



where 



A = f (C-6) 
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